% this is  Si rib waveguide 
% modes when width or thickness of core is veried
% always looking for 6 modes, incorrect solution throgh away
% can be improved by using data of previous thickness


clear all

CoreThickness=0.45:0.025:1;


    

import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
model.name('Si_rib_waveguide mode analysis');

disp('started')
%%%%%%%%%%%%%% geometry
lambda=1.55; %   wavelength in um
model.param.set('lambda', [num2str(lambda) '[um]'],'wavelength');

model.param.set('WireThickness', '0.25[um]','thickness of Si');
model.param.set('WireWidth','0.5[um]', 'width of main wire');
%model.param.set('WireWidth','0.5[um]', 'width of main wire');

model.param.set('CalculationWidth','3.8[um]', 'width of main wire');
model.param.set('AirThickness','1[um]', 'width of main wire');
model.param.set('SubThickness','1.6[um]', 'width of main wire');

geom1 = model.geom.create('geom1', 2);  
model.geom('geom1').lengthUnit([native2unicode(hex2dec('00b5'), 'Cp1252') 'm']);  % set length unit in um

% SiO2 sub
r1=model.geom('geom1').feature.create('r1', 'Rectangle');
r1.set('size',{'CalculationWidth' 'SubThickness'});
r1.set('pos',{'-CalculationWidth/2' '-SubThickness'});

% Si core
r2=model.geom('geom1').feature.create('r2', 'Rectangle');
r2.set('size',{'WireWidth' 'WireThickness'});
r2.set('pos',{'-WireWidth/2' '0'});

% air
r3=model.geom('geom1').feature.create('r3', 'Rectangle');
r3.set('size',{'CalculationWidth' 'AirThickness'});
r3.set('pos',{'-CalculationWidth/2' '0'});

geom1.run;
model.geom('geom1').feature.create('fil1', 'Fillet');
model.geom('geom1').feature('fil1').selection('point').set('r2', [3 4]);
model.geom('geom1').feature('fil1').set('radius', '0.03');

%geom1.run;

%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% materials
%n = sqrt( 1 + 0.6961663*lambda^2/(lambda^2-0.0684043^2) + 0.4079426*lambda^2/(lambda^2-0.1162414^2) + 0.8974794*lambda^2/(lambda^2-9.896161^2));

%SiO2   substrate
model.material.create('mat1');
model.material('mat1').name('SiO2');
model.material('mat1').propertyGroup.create('RefractiveIndex', 'Refractive index');
nSubstrate=nSiO2(lambda);
model.material('mat1').propertyGroup('RefractiveIndex').set('n', {num2str(nSubstrate)});
model.material('mat1').selection.set([1]);

%Si
model.material.create('mat2');
model.material('mat2').name('Si');
model.material('mat2').propertyGroup.create('RefractiveIndex', 'Refractive index');
nCore=nSi(lambda);
model.material('mat2').propertyGroup('RefractiveIndex').set('n', {num2str(nCore)});

model.material('mat2').selection.set([3]);

%air
model.material.create('mat3');
model.material('mat3').name('air');
model.material('mat3').propertyGroup.create('RefractiveIndex', 'Refractive index');
model.material('mat3').propertyGroup('RefractiveIndex').set('n', {'1'});
model.material('mat3').selection.set([2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%
% create mesh
%  create mesh
mesh1=model.mesh.create('mesh1', 'geom1');
meshSize=mesh1.feature('size');
meshSize.set('hauto', '1'); % 1-extremily fine 2 extra fine  5 normal
mesh1.feature.create('ftri1', 'FreeTri');
%mesh1.run;
%mesh1.data.transferMesh;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.physics.create('emw', 'ElectromagneticWaves', 'geom1');

model.study.create('std1');
model.study('std1').feature.create('mode', 'ModeAnalysis');

model.physics('emw').feature('wee1').set('DisplacementFieldModel', 1, 'RefractiveIndex');
model.physics('emw').feature('wee1').set('n_mat', 1, 'from_mat');

model.study('std1').feature('mode').set('geomselection', 'geom1');
model.study('std1').feature('mode').set('physselection', 'emw');
model.study('std1').feature('mode').set('shift', num2str((nSubstrate+nCore)/2));  %  aproximate mode refractive index
model.study('std1').feature('mode').set('modeFreq', ['c_const/' num2str(lambda) '[um]']);


model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'mode');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('e1', 'Eigenvalue');
model.sol('sol1').feature('e1').set('shift',  num2str((nSubstrate+nCore)/2));   %  aproximate mode refractive index
model.sol('sol1').feature('e1').set('neigs', 6);                            %  number of mode to search
model.sol('sol1').feature('e1').set('transform', 'effective_mode_index');
model.sol('sol1').feature('e1').set('control', 'mode');
model.sol('sol1').feature('e1').feature('aDef').set('complexfun', true);
model.sol('sol1').attach('std1');



for jT=1:length(CoreThickness)
     model.param.set('WireWidth',[num2str(CoreThickness(jT)) '[um]'], 'width of main wire');
     model.sol('sol1').runAll;
disp(['thickness ' num2str(CoreThickness(jT))])
clear n_eff_all mode TEmod TMmod 
% find which solutions are realistic
n_eff_all = mphglobal(model,'emw.neff','Complexout','on');
jMode=0;
for j=1:length(n_eff_all)
    if real(n_eff_all(j))>nSubstrate       
        jMode=jMode+1;
        mode(jMode).neff=n_eff_all(j);
        mode(jMode).solutionN=j;
    end
end
disp([' mode total number:  ' num2str(jMode)])

%%%%%% find whether the mode TE or TM
jTM=0;jTE=0;jTEM=0;

for j=1:jMode
    
mode(j).TEintens=mphint(model,'abs(emw.Ex)^2','Solnum',num2str(mode(j).solutionN));  % integration over all domains
mode(j).TMintens=mphint(model,'abs(emw.Ey)^2','Solnum',num2str(mode(j).solutionN));
if mode(j).TEintens>1*mode(j).TMintens 
    mode(j).Polar='TE';
    jTE=jTE+1;TEmod(jTE,1)=mode(j).solutionN;TEmod(jTE,2)=mode(j).neff;
    
elseif mode(j).TMintens>1*mode(j).TEintens
    mode(j).Polar='TM';
    jTM=jTM+1;TMmod(jTM,1)=mode(j).solutionN;TMmod(jTM,2)=mode(j).neff;
else
    mode(j).Polar='TEM';
    jTEM=jTEM+1;TEMmod(jTEM,1)=mode(j).solutionN;TEMmod(jTEM,2)=mode(j).neff;
end

end

% put each mode in order, mode with highest refractive index first
mes='';
if jTE>0 TEmod=sortrows(TEmod,-2); mes=['number of TE modes ' num2str(size(TEmod,1))];end
if jTM>0 TMmod=sortrows(TMmod,-2);mes=[mes ' of TM modes ' num2str(size(TMmod,1))]; end
 if jTEM>0 TEMmod=sortrows(TEMmod,-2);mes=[mes ' of TEM modes ' num2str(size(TEMmod,1))]; end

disp(mes )

% TE(jT,jmode)=effective refractive index  j=number corespon. thickness, j=2  effective refractive index

%%%%%%%%%%%% plot result
jplot=3;
  
if jTE>0
for jmode=1:size(TEmod,1)
   TE(jT,jmode)=TEmod(jmode,2); 
end 
end

if jTM>0
for jmode=1:size(TMmod,1)
   TM(jT,jmode)=TMmod(jmode,2); 
end 
end

end





% remove all zeros
jTE=size(TE);
for j1=1:jTE(1)
    for j2=1:jTE(2)
        if TE(j1,j2)==0;TE(j1,j2)=NaN;end
    end
end

jTM=size(TM);
for j1=1:jTM(1)
    for j2=1:jTM(2)
        if TM(j1,j2)==0;TM(j1,j2)=NaN;end
    end
end
            
% plot
figure(1); hold on
for jMode=1:jTE(2)
plot(CoreThickness, TE(:,jMode),'-r');
end

for jMode=1:jTM(2)
plot(CoreThickness, TM(:,jMode),'-b');
end
hold off

disp('OK')